%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   变量定义  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 近地点距离
h_low = 200000;

% 远地点距离
h_high = 51000000;

% 地球半径
r = 6378000;

% 轨道周期
period = 16 * 60 * 60;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   椭圆 长短半轴、焦距计算 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

a = (h_low + h_high + 2 * r) / 2;
c = (h_high - h_low) / 2;
b = sqrt(a^2 - c^2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   椭圆 轨道长度计算 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 匿名函数
f = @(x) 4 * sqrt(a^2 * cos(x).^2 + b^2 * sin(x).^2);
L = integral(f, 0, pi/2);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   近地点、平均速度计算 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
v = L / period;
area = (pi * a * b) / period;
v_max = (2 * area) / (h_low + r);

fprintf("轨道长度为： %f Km\n", L/1000)
fprintf("平均速度为： %f Km/s\n", v/1000)
fprintf("近地点速度为: %f Km/s\n", v_max/1000)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%   轨道可视化
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = linspace(0, 2*pi, 1000);
x = a * cos(t);
y = b * sin(t);

% 以30度角度旋转椭圆
angle = 30; % 角度
angle = angle / 180 * pi; % 转换后的角度
p = [x' y'];
alpha = [cos(angle) -sin(angle)
       sin(angle) cos(angle)];
p_new = p * alpha; % 坐标转换

% 动态绘图
figure(1);
x = [p_new(1, 1)];
y = [p_new(1, 2)];
m = plot(x, y, 'g');
title('卫星椭圆轨道图');
pic_num = 1;
for t = 2: 1000
    x = [x p_new(t,1)];
    y = [y p_new(t,2)];
    set(m, 'XData', x, 'YData', y)
    drawnow
    % 以下9行代码用于存储gif动图
    F=getframe(1);
    I=frame2im(F);
    [I,map]=rgb2ind(I,256);
    if pic_num == 1
        imwrite(I,map,'卫星椭圆轨道图.gif','gif', 'Loopcount',inf,'DelayTime',0.001);
    else
        imwrite(I,map,'卫星椭圆轨道图.gif','gif','WriteMode','append','DelayTime',0.001);
    end
    pic_num = pic_num + 1;
end
% saveas(1, '卫星椭圆轨道图.jpg');
